
clear;
close all;

%% 定义信道参数
%--------------------------------------------------------------------------

M_mod = 4;                         % 4QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
fc = 12500;
c = 1500;                                                     % speed of light (m/s)
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fm = dftmtx(M);       % Generate the DFT matrix
Fm = Fm./norm(Fm);    % normalize the DFT matrix

Wn=fwht(eye(N));      % Generate the WHT matrix
Wn=Wn./norm(Wn);      % normalize the WHT matrix

Md = M-Mzp;

%% select channel
P="D:\APPdocu\python\UWA_GAN\data";
GAN_Path=P;
% flag=1; % 1: NOF1; 2: NSC1
temp=1;

channel2='Generated_data\NCS1_256'; 
filepath2 = fullfile(GAN_Path, channel2, 'mat');
filename2 = fullfile(filepath2, ['mat_' sprintf('%03d', temp) '.mat']);
var_struct1=load(filename2);
name_cell=fieldnames(var_struct1);
mat2=getfield(var_struct1,char(name_cell));
h2_real=mat2(:,:,1);
h2_img=mat2(:,:,2);
h2=h2_real+1i*h2_img;
[~,G2,maxDelay2]=channel_acc(h2,M,N,2);
c2_AFDM2 = 1/(2*M);                              % Sufficiently small irrational number

%% 发送信号
    data_grid = zeros(M,N);
    data_grid(1:Md,:) = 1;
    N_syms_perfram = sum(sum(data_grid));
    
    delta_f = scs;
    T = 1/delta_f;
                 
    data_na= randi([0,1],N_syms_perfram*k,1);
    data=qammod(reshape(data_na,k,Md*N), M_mod,'gray','InputType','bit');
    X = Gen_2D_data_grid(N,M,data,data_grid);

    

%% detect
    eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
    SNR_dB = -5:5:20;
    SNR = 10.^(SNR_dB/10);
    noise_var_sqrt = sqrt(1./SNR);
    sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;     

    errorRate1 = comm.ErrorRate('ResetInputPort',true);
    berOTFS1 = zeros(length(SNR),3);
    berOTFS2 = zeros(length(SNR),3);

    errorRate2 = comm.ErrorRate('ResetInputPort',true);
    berOFDM1 = zeros(length(SNR),3);
    berOFDM2 = zeros(length(SNR),3);

    errorRate3 = comm.ErrorRate('ResetInputPort',true);
    berOTSM1 = zeros(length(SNR),3);
    berOTSM2 = zeros(length(SNR),3);

    errorRate4 = comm.ErrorRate('ResetInputPort',true);
    berOCDM1 = zeros(length(SNR),3);
    berOCDM2 = zeros(length(SNR),3);

    errorRate5 = comm.ErrorRate('ResetInputPort',true);
    berAFDM1 = zeros(length(SNR),3);
    berAFDM2 = zeros(length(SNR),3);
 


    maxDoppler1=10:10:50;
    for m = 1:length(SNR)
        disp(m)
        
        c1_AFDM1 = (2*(maxDoppler1) + 1)/(2*M); % Satisfying the orthogonality condition
        

        X1_tilda=AFDMmod(X,c1_AFDM1(1),c2_AFDM2);
        s1= reshape(X1_tilda,N*M,1);
    
        X2_tilda=AFDMmod(X,c1_AFDM1(2),c2_AFDM2);
        s2= reshape(X2_tilda,N*M,1);
    
        X3_tilda=AFDMmod(X,c1_AFDM1(3),c2_AFDM2);
        s3= reshape(X3_tilda,N*M,1);
    
        X4_tilda=AFDMmod(X,c1_AFDM1(4),c2_AFDM2);
        s4= reshape(X4_tilda,N*M,1);
    
        X5_tilda=AFDMmod(X,c1_AFDM1(5),c2_AFDM2);
        s5= reshape(X5_tilda,N*M,1);
        
         %% channel output%%%%%
    
        r1=G2*s1;
        r2=G2*s2;
        r3=G2*s3;
        r4=G2*s4;
        r5=G2*s5;
        noise= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));
        r1_n=r1+noise;
        r2_n=r2+noise;
        r3_n=r3+noise;
        r4_n=r4+noise;
        r5_n=r5+noise;

        % r12_n=r12+noise;
        % r22_n=r22+noise;
        % r32_n=r32+noise;
        % r42_n=r42+noise;
        % r52_n=r52+noise;


        [est_LMMSE_r1_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r1_n,G2,c1_AFDM1(1),c2_AFDM2); 
        [est_LMMSE_r2_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r2_n,G2,c1_AFDM1(2),c2_AFDM2);  
        [est_LMMSE_r3_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r3_n,G2,c1_AFDM1(3),c2_AFDM2);  
        [est_LMMSE_r4_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r4_n,G2,c1_AFDM1(4),c2_AFDM2);  
        [est_LMMSE_r5_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r5_n,G2,c1_AFDM1(5),c2_AFDM2); 

        % [est_LMMSE_r12_n,~] = LMMSE_detector_OTFS(N,M,M_mod,sigma_2(m),data_grid,r12_n,G2); 
        % [est_LMMSE_r22_n,~] = LMMSE_detector_OFDM(N,M,M_mod,sigma_2(m),data_grid,r22_n,G2); 
        % [est_LMMSE_r32_n,~] = LMMSE_detector_OTSM(N,M,M_mod,sigma_2(m),data_grid,r32_n,G2); 
        % [est_LMMSE_r42_n,~] = LMMSE_detector_OCDM(N,M,M_mod,sigma_2(m),data_grid,r42_n,G2); 
        % [est_LMMSE_r52_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r52_n,G2,c1_AFDM2,c2_AFDM2); 

        berOTFS1(m,:)=errorRate1(data_na,est_LMMSE_r1_n,1).';
        berOFDM1(m,:)=errorRate2(data_na,est_LMMSE_r2_n,1).';
        berOTSM1(m,:)=errorRate3(data_na,est_LMMSE_r3_n,1).';
        berOCDM1(m,:)=errorRate4(data_na,est_LMMSE_r4_n,1).';
        berAFDM1(m,:)=errorRate5(data_na,est_LMMSE_r5_n,1).';

        % berOTFS2(m,:)=errorRate1(data_na,est_LMMSE_r12_n,1).';
        % berOFDM2(m,:)=errorRate2(data_na,est_LMMSE_r22_n,1).';
        % berOTSM2(m,:)=errorRate3(data_na,est_LMMSE_r32_n,1).';
        % berOCDM2(m,:)=errorRate4(data_na,est_LMMSE_r42_n,1).';
        % berAFDM2(m,:)=errorRate5(data_na,est_LMMSE_r52_n,1).';
        
    end

%% figure
% berOTFS1(~berOTFS1)=1e-6;
% berOFDM1(~berOFDM1)=1e-6;
% berOTSM1(~berOTSM1)=1e-6;
% berOCDM1(~berOCDM1)=1e-6;
% berAFDM1(~berAFDM1)=1e-6;
% 
% berOTFS2(~berOTFS1)=1e-6;
% berOFDM2(~berOFDM1)=1e-6;
% berOTSM2(~berOTSM1)=1e-6;
% berOCDM2(~berOCDM1)=1e-6;
% berAFDM2(~berAFDM1)=1e-6;

figure

semilogy(SNR_dB,berOTFS1(:,1),'-r');             %Plot simulated BER w/ OTFS
hold on;
semilogy(SNR_dB,berOFDM1(:,1),'-b');
semilogy(SNR_dB,berOTSM1(:,1),'-g');
semilogy(SNR_dB,berOCDM1(:,1),'-x');
semilogy(SNR_dB,berAFDM1(:,1),'-o');

% semilogy(SNR_dB,berOTFS2(:,1),'-ob');  
% semilogy(SNR_dB,berOFDM2(:,1),'-+b');
% semilogy(SNR_dB,berOTSM2(:,1),'-*b');
% semilogy(SNR_dB,berOCDM2(:,1),'-.b');
% semilogy(SNR_dB,berAFDM2(:,1),'-xb');

ylabel('BER');
xlabel('SNR/dB');
% legend( 'OTFS-NOF1',"OFDM-NOF1","OTSM-NOF1","OCDM-NOF1","AFDM-NOF1",'OTFS-NSC1',"OFDM-NSC1","OTSM-NSC1","OCDM-NSC1","AFDM-NSC1");
legend( '1','11','21','31','41');
title("输出误码率")
xlim([-6,15])
grid on;
hold off;
